#####
# Lincoln Nebraska Replication
#####

rm(list=ls())
library(tidyverse)
library(broom)
library(haven)

# read in data
df <- read_csv('replication_data_lincoln.csv')

# formula
form <- gender_pct_yes ~ female_male_ratio + pct_college + pct_fem_unem + pct_minority + pct_over_65 + pct_dem + pop_density + geo_size

#run models
lm1 <- lm(form, df) 
lm1a <- lm(form, subset(df,pct_female_over_18 > median(pct_female_over_18, na.rm=T)))
lm1b <- lm(form, subset(df,pct_female_over_18 <= median(pct_female_over_18, na.rm=T)))
lm2 <- lm(pct_bond_yes ~ female_male_ratio + pct_college + pct_fem_unem + pct_minority + pct_over_65 + pct_dem + pop_density + geo_size, df) 
lm3 <- lm(pct_citizen_yes ~ female_male_ratio + pct_college + pct_fem_unem + pct_minority + pct_over_65 + pct_dem + pop_density + geo_size, df) 
lm4 <- lm(pct_coal_yes ~ female_male_ratio + pct_college + pct_fem_unem + pct_minority + pct_over_65 + pct_dem + pop_density + geo_size, df) 
lm5 <- lm(pct_audi_yes ~ female_male_ratio + pct_college + pct_fem_unem + pct_minority + pct_over_65 + pct_dem + pop_density + geo_size, df) 

# write table 
stargazer::stargazer(lm1, lm1a, lm1b, lm2, lm3, lm4,lm5,
          type='html',
          covariate.labels = c(
            'F:M Income Ratio',
            'Education Rates',
            'Pct Female Unemp',
            'Pct Minority',
            'Percent > 65yrs',
            'Percent Democrat',
            'Population Density',
            'Geographic Size'),
          star.cutoffs = c(0.10,0.05,0.01),
          dep.var.labels = c(
            'Gender', 'Bond','Draft','Coal','Auditorium'
          ), out = 'reg_table.html')

# plot results
broom::tidy(lm1) %>%
  filter(term != '(Intercept)') %>%
  mutate(term = case_when(
    term == 'female_male_ratio' ~ 'F:M Income Ratio',
    term == 'pct_college' ~ 'Education Rates',
    term == 'pct_fem_unem' ~ 'Female Unemployment',
    term == 'pct_minority' ~ 'Percent Minority',
    term == 'pct_over_65' ~ 'Percent > 65yrs',
    term == 'pct_dem' ~ 'Percent Democrat',
    term == 'pop_density' ~ 'Population Density',
    term == 'geo_size' ~ 'Geographic Size'
  ), term = factor(term, levels=c(
    'Geographic Size', 'Population Density','Percent Democrat','Percent > 65yrs',
    'Percent Minority','Female Unemployment','Education Rates','F:M Income Ratio'
  ))) %>%
  mutate(l = estimate - qnorm(0.9) * std.error,
         u = estimate + qnorm(0.9) * std.error) %>%
  ggplot(aes(term, estimate)) +
  geom_hline(yintercept=0, color='red',linetype=2) +
  geom_errorbar(aes(ymin=l, max=u),width=.1) +
  geom_point(shape=21, fill='white',size=3) +
  theme_bw() + 
  coord_flip() +
  labs(x='',y='Point Estimate') +
  ggsave(height=5, width=5,'regression_results.png')





 
      